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The formation of velocity vortices and density clusters is an intriguing phenomenon of freely 
cooling granular flows. In this work, the critical length scale Lc for the onset of instability is deter¬ 
mined via stability analysis of the linearized Navier-Stokes hydrodynamic equations of d-dimensional 
granular binary mixtures at moderate densities. In contrast to previous attempts, the analysis is 
not restricted to nearly elastic systems since it takes into account the nonlinear dependence of the 
transport coefficients and the cooling rate on the collisional dissipation. As expected from pre¬ 
vious results obtained in the very dilute regime, linear stability shows d — 1 transversal (shear) 
modes and a longitudinal (“heat”) mode to be unstable with respect to long enough wavelength 
excitations. The theoretical predictions also show that the origin of the instability is driven by the 
transversal component of the velocity field that becomes unstable when the system length L > Lc- 
An explicit expression of Lc is obtained in terms of the masses and diameters of the mixture, the 
composition, the volume fraction and the coefficients of restitution. Previous results derived in the 
limit of both mechanically equivalent particles and low-density mixtures are consistently recovered. 
Finally, a comparison with previous theoretical works which neglect the influence of dissipation on 
the transport coefhcients shows quantitative discrepancies for strong dissipation. 


I. INTRODUCTION 

Granular media are involved in many industrial and 
natural phenomena. In fact, it has been estimated that 
granular media is the second most used type of mate¬ 
rial in industry after water [l| . This is perhaps the main 
reason for which the study of granular matter has at¬ 
tracted the attention of physicists and engineers in the 
past few years. Although granular media form an ex¬ 
tremely vast family constituted by grains of different sizes 
and shapes, all these systems share relevant features. In 
particular, when granular materials are externally ex¬ 
cited (rapid flow conditions), they behave like a fluid. 
In this regime, binary collisions prevail and hence, ki¬ 
netic theory may be considered as a quite useful tool 
to describe the kinetics and hydrodynamics of the sys¬ 
tem. The main difference with respect to ordinary or 
molecular fluids is that granular systems are constituted 
by macroscopic grains that collide inelastically so that 
the total energy decreases with time. In this context, 
a granular fluid can be considered as a complex system 
that inherently is in a non-equilibrium state. In the case 
that the system is heated by an external driving force 
that compensates for the energy dissipated by collisions, 
a non-equilibrium steady state is achieved. In these con¬ 
ditions, some attempts have been recently made to for¬ 
mulate a fluctuation-response theorem based on the in¬ 
troduction of an effective temperature The gener¬ 

alization of the equilibrium fluctuation-response theorem 
to non-equilibrium states has been also confirmed in real 
experiments of intruders in driven granular fluids 0 - 
Another interesting experiments in granular matter have 
studied the response of a sheared granular medium in a 
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Couette geometry m and the behaviour of a freely ro¬ 
tating asymmetric probe immersed in a vibrated granular 
media jl2| . 

On the other hand, although significant progresses 
have been made in the past on the understanding of gran¬ 
ular flows, there are still important open challenges in 
the research of granular gases. One of the main reasons 
for which the theoretical description of these systems is 
quite intricate is that the number of relevant parame¬ 
ters needed to describe them is relatively large. This 
gives rise to a wide array of complexities that arise dur¬ 
ing the derivation of kinetic theory models. Thus, in 
order to gain some insight into the description of these 
systems in real conditions, one usually models a granu¬ 
lar fluid as a system composed by smooth hard spheres or 
disks with inelastic collisions. In this simplest model, the 
inelasticity of collisions is only accounted for by a (posi¬ 
tive) constant coefficient of normal restitution a < 1 that 
only affects the translational degrees of freedom of grains. 
Nevertheless, in spite of the simplicity of the model, it 
has been shown as a reliable prototype to explain some 
of the physical mechanisms involved in granular flows, 
especially those directly related to collisional dissipation. 

One of the most characteristic features of granular 
fluids is the spontaneous formation of velocity vortices 
and density clusters in freely cooling flows (homogeneous 
cooling state, HGS). The origin of this kind of instabil¬ 
ity is associated with the dissipative nature of collisions 
and is likely the most characteristic feature that makes 
granular flows so distinct from ordinary (elastic) fluids. 
Detected first by Goldhirsch and Zanetti [T^ and Mc¬ 
Namara [l^ in computer simulations, the instabilities 
in a free granular fluid can be well described by a lin¬ 
ear stability analysis of the Navier-Stokes hydrodynamic 
equations. This analysis provides a critical length Lc so 
that the system becomes unstable when its linear size is 
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larger than Lc- In the case of a monodisperse low-density 
granular gas, the dependence of on the coefficient of 
restitution obtained from the (inelastic) Boltzmann ki¬ 
netic equation compares quite well with numer¬ 

ical results obtained by using the direct simulation 
Monte Carlo (DSMC) method [1^. For higher densities, 
theoretical results for based on the (inelastic) Enskog 
equation [l9j | shows an excellent agreement with molec¬ 
ular dynamics (MD) simulations for a granular fluid at 
moderate density B • The stability analysis reported 
in Ref. Il9ll extends to finite dissipation some previous at¬ 
tempts |22l . [ 2 ^ carried out in the context of the Enskog 
kinetic theory but neglecting any dependence of the pres¬ 
sure and the transport coefficients on inelasticity. Nev¬ 
ertheless, while the study of the stability of the HCS has 
been widely covered in the case of granular fluids, much 
less has been made in the important subject of granular 
mixtures (namely, systems composed by grains of differ¬ 
ent masses, diameters, composition). 

Needless to say, the analysis of the stability of the HCS 
for polydiperse granular systems is much more compli¬ 
cated than for a single granular gas. Not only the num¬ 
ber of transport coefficients involved in the determination 
of the critical size Lc is higher than for a monodisperse 
gas but also they depend on more parameters, such as 
the set of coefficients of restitution characterizing the bi¬ 
nary collisions between different species. Many of the 
early attempts (23 - [?7j to obtain the Navier-Stokes coef¬ 
ficients of granular mixtures were performed by assuming 
the equipartition of granular energy. However, given that 
the lack of energy equipartition 12811 has been widely con¬ 
firmed by computer simulations (29l - l33l | and observed in 
real experiments of agitated mixtures [13, [13, the hy¬ 
pothesis of energy equipartition can only be acceptable 
for nearly elastic systems. In fact, in those previous works 
[23 - l?7j| the forms of the transport coefficients are the 
same as those obtained for ordinary mixtures [s^] and the 
influence of inelasticity is only considered in the presence 
of a sink term in the energy balance equation. A more 
rigorous derivation of linear transport f or g ranular mix¬ 
tures has been made by Garzo and Dufty [37j[ in the dilute 
regime and more recently by Garzo, Dufty and Hrenya 
|38l - l40l | for moderate densities. In these works, given that 
nonequipartition effects on transport have been consid¬ 
ered, the corresponding Navier-Stokes transport coeffi¬ 
cients exhibit an intricate nonlinear dependence on the 
coefficients of restitution of the mixture. As for single 
granular gases, the theoretical results (which have been 
obtained in the so-called fist Sonine approximation) com¬ 
pare in general quite well with computer simulations [ill- 
143 for conditions of practical interest, such as strong in¬ 
elasticity. 

The knowledge of the Navier-Stokes transport coeffi¬ 
cients of granular mixtures opens the possibility of ob¬ 
taining the critical length Lc from the (linear) stability 
analysis of the hydrodynamic equations. For dilute sys- 
tems, the theoretical predictions of kinetic theory [4^.143 
for Lc has been shown to agree very well with the DSMC 


simulations of the Boltzmann equation [ 43 ] . On the other 
hand, in spite of the explicit knowledge of the Enskog 
transport coefficients for a granular mixture [38l - l40l| . I 
am not aware of any previous solution of the linearized 
hydrodynamic equations for moderately dense granular 
mixtures. The goal of this paper is to perform a lin¬ 
ear stability analysis around the HCS in order to iden¬ 
tify the conditions for stability as functions of the wave 
vector, the volume fraction, the dimensionality of the 
system d and the parameters of the mixture (masses, 
sizes, composition and the coefficients of restitution). As 
expected, the stability analysis shows d — 1 transversal 
(shear) modes and a longitudinal heat mode to be un¬ 
stable with respect to long wavelength excitations. In 
addition, the results also show that the origin of the in¬ 
stability lies in the transversal shear mode (except for 
quite large dissipation) and hence, for sizes of the system 
larger than the critical length Lc the transversal velocity 
becomes unstable. As for dilute mixtures [13 , theoretical 
predictions for Lc compme well with recent MD simula¬ 
tions of hard spheres [49l | . A preliminary short report of 
some of the results presented here has been given in Ref. 

S. 


The plan of the paper is as follows. First, in Sec- 
tion|H]the hydrodynamic equations and associated fluxes 
to Navier-Stokes order are recalled. The explicit depen¬ 
dence of some of the transport coefficients on dissipation 
is illustrated for different systems showing that the influ¬ 
ence of inelasticity on transport is in general quite signifi¬ 
cant. Section lHll is devoted to the linear stability analysis 
around the HCS. This Section presents the main results 
of the paper. The dependence of the critical size Lc on 
the parameter space is widely investigated in section ITVl 
by varying the parameters of the system in the case of a 
common coefficient of restitution (on = 022 = 0:12 = ex). 
The paper is closed in section IVl with a brief discussion 
of the results. 


II. HYDRODYNAMIC DESCRIPTION 

We consider a binary mixture of inelastic, smooth, 
hard spheres {d = 3) or disks (d = 2) of masses mi 
and 7712, and diameters cti and (T 2 . The inelasticity of 
collisions among all pairs is characterized by three inde¬ 
pendent constant coefficients of normal restitution an, 
022 , and 012 = 021 , where is the coefficient of restitu¬ 
tion for collisions between particles of species i and j. At 
a kinetic level, all the relevant information on the state 
of the mixture is given through the one-particle velocity 
distribution function of each species fi{r,v,t) (i = 1,2). 
This quantity gives the average number of particles of 
species i that at instant t are located around the point 
r with a velocity about v. At moderate densities, the 
distribution functions fi are accurately described by the 
coupled set of inelastic Enskog kinetic equations [^ [^ . 
From this set one can derive the (macroscopic) hydrody¬ 
namic equations for the particle number density of each 
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species, 

TT-i (r,t) = y dv/,(r,v,t), (1) 

the mean flow velocity 

1 r 

U(r,t)= ^ / dvm^v/,(r,v,t), (2) 

and the granular temperature 
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r (r, t) = Y 1 J (v - U (r, t)f /,(r, v, t). 


The hydrodynamic equations are given by 


DtUi + mV ■ U + 


V • j. 


= 0, ii=l,2 


(3) 

(4) 
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FIG. 1. Plot of the (reduced) diffusion coefficient D{a)/D{1) 
as a function of the (common) coefficient of restitution a for a 
binary mixture of inelastic hard spheres (d = 3) with ui = (T 2 , 
xi = 0.5, (j> = O.l and two different values of the mass ratio 
mi/m2. 


DtV + p-^V-P = 0, (5) 


DtT--Y + A (V . q + P : VU) = -CT . (6) 
n rrii an 


In the above equations, Dt = dt +VJ ■ V is the material 
derivative, p = mini + m 2 n 2 is the total mass density, 
n = ni + n 2 is the total number density, ji is the mass 
flux for species i, q is the heat flux, P is the pressure 
tensor, and C is the cooling rate. In addition, the mass 
fluxes ji and j 2 are not independent since j 2 = —ji. 

For the two component mixture considered here there 
are d+3 independent fields, ni, n 2 , T, and U. To obtain a 
closed set of hydrodynamic equations, expressions for j^, 
q, P, and C must be given in terms of these fields. Such 
expressions are called “constitutive equations”. These 
equations have been obtained up to the Navier-Stokes 
order from the Enskog equation in Ref. [s^. They are 
given by 


m\ni 


TOiTO 2’^2 


ji =-^i^nVlnrq- 

P P 


Di 2 V In n 2 —pD^V In T, 

(7) 


q = —iVlnrii — 0^ 2 ^ \nn 2 — AVT, (8) 


Pm = pSki — T] ( ViUk + VkUi — '^^mV ■ Uj —nSkiV-V!, 

(9) 

( 10 ) 


C = + CuV • u. 


Here, Dij are the mutual diffusion coefficients, is the 
thermal diffusion coefficient, Dq^ij are the Dufour coef¬ 
ficients, A is the thermal conductivity coefficient, p is 
the hydrostatic pressure, 7] is the shear viscosity coeffi¬ 
cient, K is the bulk viscosity coefficient, and and Cu 
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FIG. 2. Plot of the (reduced) density diffusion coefficient 
Dn{a)/Dn{l) as a function of the (common) coefficient of 
restitution a for a binary mixture of inelastic hard spheres 
(d = 3) with ai — a 2 , xi = 0.5, 0 = 0.1 and two different 
values of the mass ratio mi/m^. 


are the zeroth- and first-order contributions to the cool¬ 
ing rate, respectively. The eight transport coefficients 
{Dij, Dq^ij, A, 7], K, } as well as the first-order contri¬ 
bution to the cooling rate verify a set of coupled linear 
integral equations which can be solved approximately by 
using the leading terms in a Sonine polynomial expan¬ 
sion. This type of solution provides explicit expressions 
for the complete set of Navier-Stokes transport coeffi¬ 
cients and the cooling rate in terms of the volume fraction 
(j), the coefficients of restitution and the parameters 
of the mixture (masses, sizes, and composition) [H, 

The solid volume fraction (j) is defined as 4> = (pi + 4)2 
where 


^d/2 


2‘i-idr (I) 




i = 1,2. 


(11) 


On the other hand, the expressions (Q and ([5]) for 
mass and heat fluxes can be defined in a variety of equiv¬ 
alent ways depending on the choice of the driving forces 
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FIG. 3. Plot of the (reduced) thermal diffusion coefficient 
Dria)/DriX) as a function of the (common) coefficient of 
restitution a for a binary mixture of inelastic hard spheres 
(d = 3) with = a 2 , xi = 0.5, tj) = 0.1 and two different 
values of the mass ratio mi/m 2 . 
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FIG. 4. Plot of the (reduced) shear viscosity coefficient 
ri{a)/ri{l) as a function of the (common) coefficient of resti¬ 
tution a for a binary mixture of inelastic hard spheres (d = 3) 
with (Ti = cr 2 , xi = 0.5, (j) = 0.1 and two different values of 
the mass ratio mi/m 2 . 


used. Here, to recover previous results [l^ derived in 
the stability analysis for a monodisperse dense gas, the 
hydrodynamic fields xi and n are considered instead of 
the partial densities ni and n 2 . Here, xi = n\/(n\ - 1 - 77 , 2 ) 
is the composition (or mole fraction) of species 1. With 
this simple change of variables, Eqs. o for 771 and n-i 
become 

+ - V-ji = 0, (12) 


Dtn + nV-'U^— —^V-ji=0. (13) 

77117712 

In terms of Vxi and Vti, the mass and heat fluxes read 

jr = ™ Z?„V77 - (14) 

P P T 


q = -T^D'Xxi 


-DqnVn - AVT, 


(15) 


where the transport coefficients D, D” , and Dqn are 
defined as 


D — fiDii — Di2, Dn — xi^Dii + X 2 D 12 , (16) 


D" 


Dq^l 

Xi 


Dq^2 

X2 


^qn — “f ^q,2j 


(17) 


where a ;2 = 1 — cci and /i = mi/m 2 is the mass ratio. 

The explicit expressions of the Navier-Stokes trans¬ 
port coefficients of a d-dimensional dense granular bi¬ 
nary mixture can be found in Ref. [40j |. In the case of 
mechanically equivalent particles (mi = m 2 , cti = fT 2 , 
and Uij = a), then = D” = 0 and Eqs. 

®, (HI and m agree with the results derived for a 
monodisperse dense gas [H, [s^. Moreover, in the low- 
density limit ((^ = 0), K = Cu = 0 and the results de¬ 
rived for a granular binary mixture at low-density are 
recovered [33, • Beyond the above two particular sit¬ 

uations, the dependence of the transport coefficients on 
the parameter space of the system is quite intricate. To 
illustrate this dependence, Figs. [T]® show the dimension¬ 
less quantities D{a)/D{l), Dn{a)/Dn{l), D'^{a)/D'^{!) 
and 'r]{a)/ri{l), respectively, as functions of the (com¬ 
mon) coefficient of restitution an = 022 = Q;i 2 = a for 
(j) = 0.1, ( 71/02 = 1, xi = 0.5 and two different values of 
the mass ratio m\/m 2 - Here, D{1), D„(l), and 

77(1) correspond to the values of these coefficients in the 
elastic limit. We observe that in general the influence 
of inelasticity on the transport coefficients is quite sig¬ 
nificant and so their functional form differs appreciably 
from their elastic form. This means that the previous 
predictions made for nearly elastic spheres [2J-|^ might 
quantitatively differ from those obtained here as the rate 
of dissipation increases. This will be confirmed later. 

When the expressions ®, uni), m, and (1151) for the 
pressure tensor, the cooling rate, the mass flux and the 
heat flux, respectively, are substituted into the exact bal¬ 
ance equations ® , ®, (fT^ . and (fmi) one gets the corre¬ 
sponding Navier-Stokes hydrodynamic equations for the 
hydrodynamic fields xi, 77, U, and T. They are given by 

D,:,i = . (^^^^DVxi + ™ 

77^77717772 \ P P 

+ , ( 18 ) 


D,n + ny-15 = y • ( ^^^^DVxi 

mim2 V P 

^ 77 ^^^^^^ n J, \ ^ 

P T J 


(19) 


pDtUi + Vip = Vfc 


77 (^iUk + 


+K,5kiV ■ U] 


( 20 ) 


















5 


z (a + 

/ mim2Ti 


2 __ T{m2 — mi) 

T + -pV • U =-V 

a mim2 


DVxi + + ^D'^VT 

\ P P T 

2 / \ 
f-V • T‘^D"S/xi + —A„Vn + AVT 


2 


T] ( ViUk + Vfc A - ■ U ) + 4 ^kV • U 


Vellk - nTQu^ ■ U. 


( 21 ) 


III. STABILITY OF THE LINEARIZED 
HYDRODYNAMIC EQUATIONS 


In contrast to ordinary fluids, the hydrodynamic equa¬ 
tions dH-dUl) admit nontrivial solutions even for spa¬ 
tially homogeneous states. This state is usually called 
homogeneous cooling state (HCS). In this case (no spa¬ 
tial gradients), Eqs. (fTSl) - (lri|) read 


dtXiYi = dtUYi = dtum = 0, (23) 


Note that consistency would require to consider up to sec¬ 
ond order in the spatial gradients in the expression m 
for the cooling rate, since this is the order of the terms 
appearing in the energy balance equation m coming 
from the mass flux, the pressure tensor and the heat flux. 
Thus, since the cooling rate ^ is a scalar, its most general 
form at this order for a granular binary dense mixture is 


C = + C«V • U + Cm + CmV2n2 + CtV^T 

+ CTT(Vr)^ + Cnirai(Vni)^ -|- Cn2n2 (Vna)' 
+CTm(VT) • (Vni) + CTm(VT) • (Vns) 

+ Cnira2(^^l) ■ (^^ 2 ) + iUj){'ViUj) 

+C2 ,.k(V.A)(V,A). (22 ) 


The first (linear) second-order terms (Cti^ and Cr) have 
been determined for a one-component dilute gas in Ref. 
[l5l| while all the set of coefficients (linear and nonlinear 
terms) have been computed for granular monodisperse 
gases of viscoelastic particles in Ref [s^. The evalua¬ 
tion of the above set of coefficients for granular mixtures 
is a quite intricate problem. To the best of my knowl¬ 
edge, no explicit results for these coefficients have been 
reported for granular binary mixtures, even in the sim¬ 
plest case of a low-density mixture {4> = 0). On the other 
hand, it has been shown for dilute gases that the contri¬ 
butions of the second-order terms to the cooling rate C 
are negligible [l^, as compared with the corresponding 
zeroth-order contribution C^*^^ (the first-order contribu¬ 
tion Cm vanishes for dilute gases). It is assumed here that 
the same holds in the dense case and so, for practical 
applications these second-order contributions can be in 
principle neglected in the Navier-Stokes hydrodynamic 
equations. In fact, the good agreement found in Ref 
between the present theoretical results (where the non¬ 
linear contributions to are not accounted for) and MD 
simulations for the onset of velocity vortices for strong 
inelasticity, finite density, and particle dissimilarity (see 
Figs. 2 and 3 of 0 ) supports the above expectation. 

The form of the Navier-Stokes hydrodynamic equa¬ 
tions is the same as for an ordinary binary 

mixture (a^- = 1 ), except for the presence of the contri¬ 
butions to the cooling rate (^o and and the dependence 
of the transport coefficients on the coefficients of restitu¬ 
tion. This dependence is clearly illustrated in Figs. 


dt + (a;iH, ^H, 7 h) 


Th =0, 


(24) 


where the subscript H denotes the homogeneous state. 
The dependence of the zeroth-order cooling rate = 
= C,P on xiH, riH, and Th can be estimated by 
taking Maxellian distributions for the distributions fi in 
the HCS. Here, is the partial cooling rate associated 
with the partial temperature T^, which is a measure of 
the mean kinetic energy of particles of species i. The 
expression of in the Maxwellian approximation is 


AO) 


dA2 


2 


^ ^oXi] 

i=i 


CT 12 ) I e,e, ) 


X (I -I- aij) 



fii + m) 


9i + 9j 


(25) 


where Xij is the pair distribution function at contact, 
Pij = rriiKm^ -|- m^), 9, = miT/mTi, Uij = {ui A (Xj)l2, 
m = (mi -I- 7712)12., and 


V = 


^{d-l )/2 8 jY 

r(|) d + 2^'^12 Y = 


(26) 


is an effective collision frequency chosen to recover pre¬ 
vious results found for a dense one-component granular 
gas 0 . Note that for mechanically equivalent particles 
(mi = m 2 , (Ti = 0 - 2 ), the collision frequency v is associ¬ 
ated with the elastic shear viscosity of a dilute gas. Due 
to inelasticity in collisions, energy equipartition is bro¬ 
ken and so, in general Ti 7 ^ T 2 . The temperature ratio 
7 = T 1 /T 2 is determined from the condition Ci°^ = 

In order to solve Eq. (IMl) it is convenient to change to 
a new time variable defined as 


T = 





(27) 


In this new time variable, the integration of Eq. (l24l) 
yields 


T{t) = T(0)e"^'^“^ (28) 

where Cq = jv. To find the relation between the “in¬ 
ternal” time (related to the average number of collisions 
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suffered per particle) and the “external” time t, one in¬ 
tegrates the relation for dr using v ~ y/T and gets the 
usual Half’s law 


Wk(T) = 


^Uk(T) 


6'k(T) = 


^Tk(r) 


(32) 


'Ch(t) ’ ’ 

where Syupir) = {pi.k('r), pk('r), Wk(T), dklr)} is defined 
as 


m 


no) 

(l+i^(o)(o)i)2’ 


(29) 


where (0) is the cooling rate at the initial time. The 
partial temperatures Ti also have the same time depen¬ 
dence (l29l) but each with a different value 

On the other hand, computer simulations [3 [3 have 
clearly shown that the HCS is unstable with respect to 
long enough wavelength perturbations. To analyze this 
problem it is convenient to carry on a (linear) stability 
analysis of the nonlinear hydrodynamic equations (1181) - 
(121) with respect to the homogeneous state for small ini¬ 
tial excitations. For ordinary fluid mixtures such per¬ 
turbations decay in time according to the hydrodynamic 
modes of diffusion (shear, thermal, mass) and damped 
sound propagation. The perturbation analysis made here 
is for fixed coefficients of restitution different from unity 
in the long wavelength limit. It will be seen that the cor¬ 
responding hydrodynamic modes for a granular mixture 
|47j | differs from those obtained for normal mixtures. In 
addition, as has been widely explained in some previous 
papers [3, El) 113) Sll , the linearization of Eqs. (IT^ - (I^ 
about the homogeneous base state leads to a set of cou¬ 
pled partial differential equations with coefficients that 
depend on time since the HCS is cooling. This time de¬ 
pendence can be eliminated through convenient changes 
in the time and space variables and a scaling of the hy¬ 
drodynamic fields. 

Let 5yp{T,t) = yp{r,t) — yiip{t) denote the deviation of 
{a;i,n,U,T} from their values in the HCS. If the initial 
spatial perturbation is sufficiently small, then for some 
initial time interval these deviations will remain small 
and the hydrodynamic equations (IT8) - (l2T) can be lin¬ 
earized with respect to 5yy{Y,t). As said before, to elim¬ 
inate the time dependence we introduce the time variable 
T (defined in Eq. (E71) ') and the space variable 


t 


1 ^'h(0 
2uh (t)*^’ 


(30) 


<5yk/3(T) = j Ote *'""^(52//3(£,t). (33) 

Note that in Eq. (1551) the wave vector k is dimensionless. 
In terms of these variables, the d — 1 transverse veloc¬ 
ity components Wkj_ = Wk — (wk • k)k (orthogonal to 
the wave vector k) decouple from the other four modes 
and hence can be obtained more easily. Their evolution 
equation is 


where 



Co + 2 ^**' 


Wk± = 0, 


* _ VYim 
' — 2 ' 
Ph^h 


(34) 


(35) 


In Eq. (l34l) it is understood that Co is also evaluated in 
the HCS. The solution to Eq. (|34l) is 


Wk_L(k, r) = Wki(O) exp[s_L(fc)r], (36) 


where 


snk)=Q-\ri*k^- (37) 

This identifies d—\ shear (transversal) modes analogous 
to the elastic ones [13 • According to Eq. (|57) . there 
exists a critical wave number given by 


k 


C _ 

_L — 



(38) 


This critical value separates two regimes: shear modes 
with k > always decay while those with k < k‘^ grow 
exponentially. 

The remaining (longitudinal) modes correspond to the 
composition field pi.k, the density field pk, the longitudi¬ 
nal component of the velocity field Wk|| = Wk-k (parallel 
to k), and the temperature field ^k. These modes are 
coupled and obey the time-dependent equation 


where vii{t) = \/TH(t)/m. According to Eq. (l27l) . the 
dimensionless time scale r is therefore an average number 
of collisions per particle in the time interval between 0 
and t. The unit length introduced in Eq. (1501) 

is proportional to the effective time-independent mean 
free path l/rinan^- 

A set of Fourier transformed dimensionless variables 
are then defined by 


felk(T~) 

a;iH 


Pk{r) 


fak('r) 

riH 


(39) 

where now Szkpir) denotes the four variables 
(pi.k, Pk, 6*k, rck||) • The matrices in Eq. (159) are 




/ 0 0 0 0 \ 

0 0 0 0 

-2^1 S -2^^ -Co 0 

Vo 0 0 Co/ 


( 40 ) 


Pi,k('r) 


(31) 
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= 


( 0 
0 
0 

(l+/j)a dy* 
\ 2{l+fiS) dxi 


0 
0 
0 

(l+^)<5 


din^p*) 


2 x 2 { 1 +P^S) dn* 


0 

0 

0 

(1+Ai)<5 ^ 

2x2(l+pS)^ 


0 

-1 


-h* - c. 


(41) 


m(2) = 
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(42) 


where the subscript H has been omitted in Eqs. (IM)) - (|42|l 
for the sake of brevity. In these equations, n* = naf 2-1 ^ = 
xi 1x2-1 and I have introduced the dimensionless transport 
coefficients 


^ _ mim2r'H 


D* = 


PuTi 


H 


Dl = 


, _ TOim2J^H 


PbT] 


Dn 


H 


(43) 


^ d"* = 


nnTn 


4 mi/fj 
d{d + 2) uh 


Zl", (44) 


d(d + 2) UH = 


mi/H 


d(d + 2) uhTh 


Ah, 

(45) 


K* = -n-KH- (46) 

PH'Cl 

In addition, the (reduced) hydrostatic pressure p* = 
p/{nT) is given by [s^ 


* 1 I ^ 

V =l + ^IT7aT« 2^2:, 


dr(f)'" ^-■^'■^'(^ 12 ) 


Mi* (1 + Oiij) 


where the temperature ratios 7 ^ = Ti/T are defined as 
7 1 


7i = 


l + a;i( 7 -l)’ l + xi( 7 -l)' 


(47) 
as 

(48) 


Note that p* depends explicitly on n* and xi and also 
through its dependence on Xij and 7 ^. 

The longitudinal four modes have the form exp[sn(fc)T] 
for n = 1 , 2 ,3,4, where s-n{k) are the eigenvalues of the 
matrix 


M = (49) 

In other words, they are the solutions of the quartic equa¬ 
tion A{k, s) = 0 where 

A(fc, s) = det(M — si). (50) 

Here, 1 denotes the identity matrix. Although the ex¬ 
plicit form of Aik, s) can be easily obtained from Eqs. 


(Uni-dill), its expression will be omitted here for the sake 
of brevity. 

It is instructive to consider first the solutions to Eqs. 
(E3 and (1501) in the extreme long wavelength limit, fc = 0. 
In this case, the eigenvalues of the hydrodynamic modes 
are given by 

s± = ^Co, s„ = (0,0,-a,Co)- (51) 

Two of the eigenvalues are positive, corresponding to 
growth of the initial perturbation in time. Thus, some 
of the solutions are unstable. The two zero eigenval¬ 
ues represent marginal stability solutions, while the neg¬ 
ative eigenvalue gives stable solutions. For general ini¬ 
tial perturbations all modes are excited. These modes 
correspond to evolution of the fluid due to uniform per¬ 
turbations of the HCS, i.e. a global change in the HCS 
parameters. The unstable modes are seen to arise from 
the initial perturbations u'kj_(0) or r(jj(.||(0). The unstable 
modes may appear trivial since they are due entirely to 
the normalization of the fluid velocity by the time depen¬ 
dent thermal velocity vnit). However, this normalization 
is required by the scaling of the entire set of equations to 
obtain time independent coefficients. 

For k ^ 0, the dependence of the longitudinal modes 
on the wave vector k is quite intricate. On the other 
hand, the critical longitudinal mode can be obtained 
from the equation A{k, s) = 0 when s = 0. This leads to 
the quartic equation 

Ao + A2k^ + Aik^ = 0, (52) 

where the coefficients Ai are known functions of the solid 
volume fraction, the coefficients of restitution, and the 
parameters of the mixture. The critical value fcj) is the 
largest real root of Eq. dSH). As before, given that its 
explicit expression is very long and not relevant for the 
purposes of this paper, I shall omit the form of fcj) for the 
sake of simplicity. 

For mechanically equivalent particles (mi = m 2 = m, 
CTi = (72 = cr, and aij = a), the results derived in this 
section agree with those obtained before for monodis- 
perse dense gases [T^ . Moreover, in the low-density limit 
(n* —>■ 0 ), one recovers the dispersion relations obtained 
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FIG. 5. The critical length scale Lc for velocity vortices in units of the diameter cr as a function of the coefficient of restitution 
a for a system of mechanically equivalent particles. Three different values of the volume fraction (j) are considered: (a) 0 = 0.1, 
(b) (j) = 0.2, and (c) rf) = 0-4. The left panel corresponds to disks {d — 2) while the right panel refers to spheres (d = 3). In 
each case, the system is linearly stable for points below the corresponding curve 


in a previous work for a dilute granular binary mix¬ 
ture. These limits show the consistency of the linear sta¬ 
bility analysis carried out in this paper for dense granular 
binary mixtures. 


IV. CRITICAL SIZE FOR THE ONSET OF 
INSTABILITIES 


In a system with periodic boundary conditions, the 
smallest allowed wave number is 27r/L, where L is the 
largest system length. Hence, for given values of density, 
coefficients of restitution and parameters of the mixture, 
we can identify a critical size Lc so that the system be¬ 
comes unstable when L > Lc- The value of Lc is deter¬ 
mined by equating 

^ = max{fc^, A;|j}, L* = ^^Lc- (53) 

On the other hand, the present results show that in gen¬ 
eral and hence, the origin of instability is as¬ 

sociated with the transversal components of the veloc¬ 
ity field. This result has been verified for a low-density 
monodisperse granular gas [13 and for a dilute granu¬ 
lar binary mixture by numerically solving the Boltz¬ 
mann equation by means of the Direct Simulation Monte 
Carlo (DSMC) method [l^. Recent MD simulations 
[13) mi Ii3| for granular fluids at moderate densities have 
also found that ^ latter case, the theo¬ 

retical predictions for the critical size compare well with 
MD simulations even for strong dissipation [i^]- Thus, 
according to Eqs. (1551) and (I55)) . the critical length scale 
for velocity vortex instability is given by 


0'12 


d + 2r(f) 

2V2 



(54) 


Equation (IMl) gives the dependence of the critical 
length on the parameter space of the system. Specifi¬ 
cally, the parameter space is the mass ratio mil m 2 , the 
composition ni/(ni -I- n 2 ), the ratio of diameters CTi/cr 2 ) 
the coefficients of restitution a^-, and the solid volume 
fraction (j). According to Eq. (l54l) . the dependence of Lc 
on the parameters of the mixture is essentially accounted 
for by the term i/rfjQ. The expression of the cooling 
rate Co is given by Eq. (1551) while the explicit form of the 
(reduced) shear viscosity rj* is provided in the Appendix 
HI for the sake of completeness. 

Given that the parameter space of the problem is large, 
in order to reduce the number of independent parame¬ 
ters the simplest case of a common coefficient of resti¬ 
tution (a = ail = OL 22 = 0 : 12 ) is considered. Thus, 
once the dimensionality of the system is fixed, the pa¬ 
rameter space is reduced to five dimensionless quantities: 
{mi/m2,<7i/a2,xi,cj),a}. 

Three different values of the overall volume fractions 
have been considered: (j) = O.I, (j) = 0.2, and (p = 0.4. The 
first two values of (p represent a granular fluid with mod¬ 
erate density while the latter one corresponds to a system 
with high density. Three different values of the common 
coefficient of restitution have been analyzed: a = 0.9 
(weak dissipation), a = 0.8 (moderate dissipation), and 
a = 0.5 (strong dissipation). 

Before considering a binary mixture, the case of me¬ 
chanically equivalent particles (monodisperse dense gas) 
is illustrated for inelastic hard disks {d = 2) and spheres 
(d = 3). Figure [5] shows Lcja versus the coefficient of 
restitution a for different values of the volume fraction 
p. First, we observe that although the dependence of the 
critical length on dissipation is quite similar in d = 2 and 
d = 3, the magnitude of Lc (measured in units of the 
diameter cr) is larger in disks than spheres. Thus, for a 
given value of a, the critical length for velocity vortices 
increases as the dimensionality of the system decreases. 
With respect to the dependence on the density p, at a 
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FIG. 6. The critical length scale Lc for velocity vortices in units of ai 2 as a function of the mass ratio mijm 2 . We have 
considered a granular binary mixture with xi = 0.5, cri/cr 2 = 1 and (/> = 0.1. Three different values of the (common) coefficient 
of restitution are studied: (a) a = 0.9, (b) a = 0.8, and (c) a = 0.5. The left panel corresponds to disks (d = 2) while the right 
panel refers to spheres (d = 3). In each case, the system is linearly stable for points below the corresponding curve. 




FIG. 7. The critical length scale Lc for velocity vortices in units of g \2 as a function of the ratio of diameters (Ti/(T 2 . We 
have considered a granular binary mixture with xi = 0.5, mi/m 2 = 5 and </> = 0.2. Three different values of the (common) 
coefficient of restitution are studied: (a) a = 0.9, (b) a = 0.8, and (c) a = 0.5. The left panel corresponds to disks (d = 2) 
while the right panel refers to spheres (d = 3). In each case, the system is linearly stable for points below the corresponding 
curve. 



FIG. 8. The critical length scale Lc for velocity vortices in units of ai 2 as a function of the concentration xi. We have considered 
a granular binary mixture with mi/m 2 = 5, cri/(T 2 = 1 and (/> = 0.2. Three different values of the (common) coefficient of 
restitution are studied: (a) a = 0.9, (b) a = 0.8, and (c) a = 0.5. The left panel corresponds to disks (d = 2) while the right 
panel refers to spheres (d = 3). In each case, the system is linearly stable for points below the corresponding curve. 
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FIG. 9. The (dimensionless) critical length scale Lc/avi as 
a function of the mass ratio m\lmi for a binary mixture of 
inelastic hard spheres (d = 3) with x\ =0.1, axja^ = I, (fi = 
0.1 and a (common) coefficient of restitution a = 0.7. The 
solid fine is the result derived here while the dashed line is the 
result derived from Eq. dSl by assuming energy equipartition 
(q- = 1 ) and by neglecting the a-dependence of 77 *. 


given value of dissipation, it is quite apparent that 
decreases with density and hence, smaller systems are 
required to observe the shearing instability as the gran¬ 
ular mixture becomes denser. 

Now we consider granular binary mixtures. Figure [5] 
shows Lc/o ’12 as a function of the mass ratio raxjrai with 
xi = 0.5, cri/o -2 = !,</’ = 0.1 and different values of a. 
Regarding the influence of dimensionality on the critical 
length, for given values of the coefficient of restitution 
and the parameters of the mixture, the critical length is 
(significantly) larger for c? = 2 than for d = 3. As will 
show later, this trend is also observed in the remaining 
plots presented in this paper. Moreover, at a given value 
of mi/m 2 , as expected Lc decreases with collisional dis¬ 
sipation. In fact, Lc —f 00 for elastic collisions (a = 1). 
The theoretical results also show that the Lc/(Ti 2 predic¬ 
tions for a = 0.8 and a = 0.5 seem to converge for large 
mass ratios and eventually crossover. This tendency is 
consistent with MD simulations (see for instance, fig. 2 
of Ref. [H). Figure [7] shows the critical size as a function 
of the ratio of diameters with Xi = 0.5, mi/m 2 = 5 and 
(j) = 0.2. It is quite apparent that Lc/<Ji 2 increases with 
the ratio (Ji/a 2 , except in a small region of ai /02 close 
to one. In addition, the influence of the size ratio on the 
critical length is stronger for weak dissipation (a = 0.9) 
than for strong dissipation (a = 0.5), this behavior being 
independent of the dimensionality of system. Finally, the 
dependence of the critical length on species composition 
xi is illustrated in Fig. ([5]) for a binary mixture with 
o’i/cr 2 = I, mi/m 2 = 5 and cj) = 0.2. We see that Lc is 
not a monotonic function of xi. Moreover, the impact 
of composition on the critical size is much less important 
than the one observed with respect to the mass ratio 
(see Fig. ini) and/or the size ratio (see Fig. [7]). In fact, for 
a = 0.8, Lc displays a very weak dependence on xi. It 
is important to recall that the theoretical predictions for 
Lc derived from the (inelastic) Enskog kinetic equation 


have been recently assessed against MD simulations of 
inelastic hard spheres [1^. The comparison carried out 
in this work has shown in general an excellent agreement 
between theory and simulation when physical properties 
of particles are not quite disparate. This agreement be¬ 
comes only good (at worst 20% error) for the most ex¬ 
treme conditions studied [see for instance. Fig. 2(b) for 
toi/to 2 = 10, xi = O.I, (f) = 0.2 and a = 0.7]. On the 
other hand, based on the results obtained for the tracer 
diffusion coefficient one would expect that the ac¬ 

curacy of the Enskog result for the shear viscosity (which 
is the transport coefficient involved in the evaluation of 
the critical length scale) by considering only the leading 
term in a Sonine polynomial expansion |39j would de¬ 
cease as the mass ratio becomes more disparate. Thus, 
it is reasonably to expect that the second Sonine correc¬ 
tion to rj mitigate part of the discrepancies observed in 
Ref. [ 4 ^ for high dissipation and/or extreme mass or size 
ratios. 

Before closing this Section, it is interesting to com¬ 
pare the predictions offered here with those obtained be¬ 
fore [23 - 1^ by neglecting the nonlinear dependence of 
the transport coefficients on dissipation. As mentioned 
in the Introduction, these works [2J-[^ are devoted to 
nearly elastic systems and hence, they assume the valid¬ 
ity of energy equipartition. In fact, according to this level 
of approximation, the inelasticity in collisions is only ac¬ 
counted for in the cooling rate Co and the expressions 
of the Navier-Stokes transport coefficients are the same 
as those obtained for elastic collisions [s^. In order to 
gauge the impact of the non-equipartition of granular en¬ 
ergy and the a-dependence of the transport coefficients. 
Fig. [5] compares the dependence of L^ on mi/m 2 given 
by the present calculations (solid fine) with those de¬ 
rived before (dashed line) (234^ . In this plot, xi = 0.1, 
ai/(T 2 = 1, C* = 0-1 and a = 0.7. The comparison in¬ 
dicates good agreement between both approaches in the 
monodisperse gas case (mi = m 2 ). On the other hand, 
quantitative significant differences appear as the mass ra¬ 
tio increases. 


V. SUMMARY AND DISCUSSION 

Particulate or granular flows play a critical role 
in chemical process industries (e.g., fluidized catalytic 
cracking), energy (e.g., gasification of coal and biomass), 
pharmaceuticals (e.g., powder processing, granulation), 
geological phenomena (e.g., planetary formation), and 
agriculture (e.g., grain conveying). In spite of this per¬ 
vasiveness, granular systems are not completely under¬ 
stood yet. Apart from its practical interest, the study of 
granular matter under fluidization conditions poses open 
challenges from a fundamental point of view. In fact, 
granular matter can be considered as a good example 
of a complex system since it is composed of many de¬ 
grees of freedom (many particles) that are imbedded in 
a network of strong nonlinear interactions. This aspect 




11 


explains in part because the use of non-equilibrium statis¬ 
tical physics or kinetic theory to describe granular fluids 
(i.e., when the material is externally excited) has been an 
active area of research in the past several decades. On 
the other hand, although in many conditions the motion 
of grains exhibits a great similarity to the random mo¬ 
tion of atoms or molecules of an ordinary gas, the fact 
that collisions between grains are inelastic gives rise to 
subtle modifications of the conventional hydrodynamic 
equations. 

In this paper a granular fluid mixture has been mod¬ 
eled as a mixture composed by smooth, inelastic hard 
spheres. The simplest situation for this system cor¬ 
responds to the so-called homogeneous cooling state 
(HCS). It describes a uniform state with vanishing flow 
field and a granular temperature decreasing monotoni- 
cally in time. However, it is well-known that the HCS 
can be instable against long wavelength spatial pertur¬ 
bations, leading to cluster and vortex formation. The 
existence of the instability can be predicted via a linear 
stability analysis of the nonlinear Navier-Stokes hydro- 
dynamic equations. This analysis allows us to identify a 
critical length size Lc beyond which the system becomes 
unstable. On the other hand, given that the expression of 
Lc involves the Navier-Stokes transport coefficients of the 
system, one has to explicitly determine first these coeffi¬ 
cients to get the impact of the parameter space (masses 
and sizes of grains, coefficients of restitution, composi¬ 
tion, density, • • •) on the critical size. 

In order to gain some insight into the general problem, 
a kinetic theory description has been adopted where all 
the relevant information on the state of the system is 
given through the knowledge of the one-body distribu¬ 
tion function for each species. For moderate densities, 
the Enskog equation [^, for smooth inelastic hard 
spheres can be considered as a reliable kinetic equation. 
As in the case of elastic collisions, the Enskog equation 
neglects velocity correlations between the particles which 
are about to collide (molecular chaos) but takes into ac¬ 
count spatial correlations. To first order in the spatial 
gradients, the set of Enskog equations of the mixture have 
been solved [ssl - l^ by means of the Chapman-Enskog 
method conveniently adapted to dissipative dynam¬ 
ics. As for ordinary fluid mixtures [s^, the Navier-Stokes 
transport coefficients are given in terms of a set of linear 
integral equations [s^ that can be approximately solved 
by considering the leading terms [s^, in a Sonine 
polynomial expansion. The knowledge of the transport 
coefficients allows one to obtain the critical size for the 
onset of instability in terms of the parameters of the mix¬ 
ture. The present paper has addressed this problem for a 
d-dimensional granular binary mixture described by the 
Enskog equation. The results derived here cover some of 
the aspects not accounted for in previous studies. Specif¬ 
ically, (i) it takes into account the nonlinear dependence 
of transport on collisional dissipation (and thus, the the¬ 
ory is expected to hold for a wide range of coefficients of 
restitution); (ii) it considers the influence of energy non- 


equipartition on the critical size; and (iii) it is applicable 
to moderate densities (say for instance, volume fraction 
typically smaller than or equal to 0.25 for hard spheres). 
Thus, the theory subsumes all previous studies for single 
fluids [l^ [l^ [12, nil , dilute mixtures [H, jil] and nearly 
elastic dense mixtures [2l - l27j| , which are recovered in the 
appropriate limits. 

Our findings agree qualitatively well with previous re¬ 
sults m, HI] carried out by using the elastic expressions 
of the Navier-Stokes transport coefficients since the ef¬ 
fect of dissipation on transport coefficients do not sig¬ 
nificantly change the form of the dispersion relations. 
Thus, stability analysis of the linearized hydrodynamic 
equations for the mixture shows d—1 transversal (shear) 
modes and one longitudinal heat mode. An analysis of 
the dependence of the above modes on the parameter 
space shows that in general the instability is driven by 
the transversal shear mode and hence, the critical size is 
given by Eq. (I5i)) where the (reduced) cooling rate Q is 
defined by Eq. (1^ while the expressions defining the (re¬ 
duced) shear viscosity rj* are provided in the Appendix. 
It is quite apparent that Ldoyi presents a complex de¬ 
pendence on the parameter space of the system so that, 
it is intricate to disentangle the impact of the different 
parameters on the critical size for instability. Thus, to re¬ 
duce the number of independent parameters, a common 
coefficient of restitution (a^- = a) has been assumed to 
illustrate the influence of the mass and size ratios as well 
as the composition on Lc- The results have been car¬ 
ried out in Section [Tv] for hard disks {d = 2) and spheres 
(d = 3) in Figs. [SHHl We observe that the role played by 
the mass and size ratios on the critical size is more rele¬ 
vant than that of the composition. In addition, although 
the dependence of Lc on the parameter space is qualita¬ 
tively similar for disks and spheres, the magnitude of the 
critical size is larger in the case of two dimensions than 
three dimensions for the same system. 

With respect to the effect of collisional dissipation on 
Lc, as expected quantitative discrepancies between the 
present results and those obtained by assuming energy 
equipa rtition and elastic forms for transport coefficients 
[241 - 1^ appear as the coefficient of restitution decreases. 
This is clearly illustrated in Fig. |9l Therefore, although 
the results derived in the quasielastic limit could predict 
reasonably well the dispersion relations, one expects that 
the results reported here improve these findings for condi¬ 
tions of practical interest where the transport coefficients 
are clearly affected by dissipation. 

Although the present theory applies in principle for ar¬ 
bitrary values of the coefficients of restitution, it has some 
important restrictions. First, given that some previous 
computer simulation works [59|-|6l| have clearly shown 
that the molecular chaos hypothesis fails for inelastic col¬ 
lisions as the density increases, it is possible that the 
limitations of the Enskog equation are greater than for 
elastic collisions. However, as mentioned in the Introduc¬ 
tion, despite the above limitation the theoretical predic¬ 
tions obtained here from the Enskog equation compare 
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well with recent MD simulations of hard spheres. 
This shows again the reliability of the Enskog theory 
to accurately describe macroscopic properties (such as 
transport coefhcients and/or the onset of instability) for 
a wide range of densities and/or coefficients of restitu¬ 
tion. Another important limitation of the theory is that 
the critical length obtained here has been estimated by 
considering only the first Sonine approximation for the 
shear viscosity coefficient 77 . Recent results [H, |4^ for 
the tracer limit have shown that the reliability of the first 
Sonine solution can be questionable for strong collisional 
dissipation and/or disparate values of the mass and diam¬ 
eter ratios. Therefore, the second Sonine correction to 77 
could in part improve the agreement between theory and 
MD simulations [i^ for Lc when the mass ratio is large. 
The evaluation of the second Sonine approximation to 
the shear viscosity is an interesting open problem to be 
studied in the near future. Another possible direction of 
study is to consider the so-called modified Sonine method 
[6^[63j to determine 77 . This new approach consists of re¬ 
placing, where appropriate in the Chapman-Enskog pro¬ 
cedure, the MaxwellBoltzmann distribution weight func¬ 
tion (used in the standard first Sonine approximation) 
by the homogeneous cooling state distribution for each 
species. As in the case of dilute binary mixtures [b^ . 
it is expected that the modihed Sonine approximation 
improves the estimates of the standard one at strong dis¬ 
sipation. Work along these lines will be carried out in 
the near future. 
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where the (reduced) bulk viscosity coefficient is 
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The partial kinetic contributions 77 /* read 
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The (reduced) collision frequencies are given by [H, 
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Appendix A: Shear viscosity of a dense granular 
binary mixture 

The explicit dependence of the (reduced) shear viscos¬ 
ity 77 * defined by Eq. (1351) on the parameter space of the 
system is provided in this Appendix. The expression for 
77 * reads 


47r(d-l)/2 
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(Al) 


The collisional contribution 77 ! to the shear viscosity is 
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The expressions for r 22 and r^i can be obtained from Eqs. 
(1X71) and (IA 8 I) . respectively, by just making the changes 
1 o 2. 

In order to get the dependence of the shear viscosity 
on the parameters of the system, one needs to know the 
explicit forms of the pair correlations functions Xij- For 
hard disks (d = 2 ), a good approximation for the pair 
correlation function Xij is 


1^9 (p o'iajMi 
^ l-p ^ 16(l-</)2 aijM2 


(A9) 
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where 


correlation function Xij the following approximation 


2 

M„ = ^ Xscr". (AlO) 

S = 1 

In the case of hard spheres (d = 3), we take for the pair 


1 3 (j) aiajM2 1 / ai(jjM2\^ 
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